Large-Scale Simulations of the Two-Dimensional Melting of Hard Disks 
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Large-scale computer simulations involving more than a million particles have been performed 
to study the melting transition in a two-dimensional hard disk fluid. The van der Waals loop 
previously observed in the pressure-density relationship of smaller simulations is shown to be an 
artifact of finite-size effects. Together with a detailed scaling analysis of the bond orientation order, 
the new results provide compelling evidence for the Halperin-Nelson- Young picture. Scaling analysis 
of the translational order also yields a lower bound for the melting density that is much higher than 
previously thought. 

PACS numbers: 64.60.Fr, 64.70.Dv 



A system of hard disks in two dimension (2D) is one of 
the simplest models of a classical fluid. But beneath the 
deceptive simplicity of this model, 2D hard disks exhibit a 
set of surprisingly rich behaviors. Unlike in three dimen- 
sions, a 2D solid possesses only quasi-long-range trans- 
lational order which decays algebraically to zero at large 
distances Q . Instead of the usual first-order transition in 
three dimensions, a 2D solid is also expected to melt into 
a liquid via two continuous transitions. The intervening 
phase called the "hexatic" was predicted by Halperin and 
Nelson 0, |j| and Young Q to possess quasi-long-range 
bond orientation order but no long-range translational 
order. 

Given the simplicity of the hard disk model, it would 
seem easy to either prove or disprove the Halperin- 
Nelson- Young (HNY) theory by detailed computer sim- 
ulation studies. But twenty- five year after the HNY the- 
ory was first proposed, simulations that could definitively 
identify the nature of the melting transition are still lack- 
ing 5] . The first simulation of 2D hard disks was carried 
out by Alder and Wainwright @. Based on the appear- 
ance of a van der Waals loop in the pressure, they con- 
cluded that the melting transition must be first-order. 
Since then, as more computing power has become avail- 
able, simulations have been carried out with increasi ngly 
larger system sizes [| BBIj P El H 13 El H US 
lialiililSlMSl2ilMl23,but instead of clarifying 
the picture, these simulations have provided conflicting 
conclusions about the nature of melting transition. One 
consensus that did emerge from the more recent simu- 
lation studies is that the 2D hard disk system is very 
sensitive to finite-size effects near the melting transition. 
This is not unexpected if the transition is continuous, 
but compared to fluids with a soft potential |2^| the hard 
disk system is much more prone to finite-size errors and 
boundary effects. In a simulation of up to N = 128 2 par- 
ticles, Zollweg and Chester observed that the equi- 
libration time increased dramatically for densities very 
close to the melting transition - systems of this size were 
apparently not large enough to reach the scaling limit. 



The largest simulation that has been performed to date 
was carried out by Jaster with up to N — 256 2 particles 
[lfl Il7j , and more recently for two higher densities with 
up to N = 1024 2 ji^l- Even though a van der Waals loop 
was observed in the pressure at densities between p = 
0.895 and 0.910 (solid squares in Fig. 1), Jaster showed 
using a scaling analysis that his data were also compatible 
with the HNY scenario. A van der Waals loop is often 
the sign of a first-order transition, but it may also arise 
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FIG. 1: Pressure of the hard disk fluid as a function of den- 
sity. Solid squares are N = 256 2 data from Jaster [ttTI . and 
crosses and plus, respectively, are N = 512 2 and N = 1024 2 
data from Jaster |2^l . Open circles and open squares are data 
from the present work for A^ = 512 2 and 1024 2 , respectively, 
with error bars as indicated. The dotted line is a guide to the 
eye through Jaster's data for N — 256 2 . The dashed and solid 
line are guides to the eye through the N — 512 2 and 1024 2 
data in the present work. Note the presence of an apparent 
van der Waals (vdW) loop between p = 0.895 and 0.910, which 
becomes shallower for increasingly larger size simulations. For 
N = 1024 2 , the van der Waals loop between p = 0.895 and 
0.905 has disappeared completely, with a small decrease in 
the pressure still visible for p — 0.910. The two arrows in- 
dicate the approximate locations of the isotropic-hexatic and 
hexatic-solid boundaries. 
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from finite-size errors. To definitively rule out a first- 
order scenario, one must demonstrate that the van der 
Waals loop is a finite-size artifact, i.e. it must be shown 
to disappear with larger size simulations. Curiously, the 
same van der Waals loop was observed for two different 
sizes in Jaster's data - the pressure for N = 128 2 (not 
shown in Fig. f ) and 256 2 coincide almost completely. 

In this letter, we describe a Monte Carlo study of 2D 
hard disks for up to N = 1024 2 = 1048576 particles. The 
calculations were carried out in the canonical ensemble, 
in a square box with periodic boundary condition and 
in a rectangular box with aspect ratio v3 : 2 for the 
higher densities. We worked with densities in the range 
p = 0.880 to 0.920, which accordingtoprevious estimates 
should span the transition region [1, la, 0, EH E3] • Den- 
sities p are given in reduced units where the hard disk 
diameter is one. 

While the rationale for going to larger system size is to 
eliminate finite-size effects, larger simulations also take 
longer to equilibrate. We have focused on N = 512 2 
to try to carry out detailed simulations covering a large 
range of densities between p = 0.880 and 0.920. At this 
size, one run at each density took several months of CPU 
time. Additional larger simulations with N — 1024 2 were 
performed for four densities between p = 0.895 and 0.910 
in the vicinity of the van der Waals loop previously ob- 
served in smaller simulations. In contrast, Jaster's recent 
simulations focuses on a different region in the phase 
diagram, offering data for p = 0.918 at N = 1024 2 and 
two densities, p = 0.914 and 0.918, for N = 512 2 . 

Two different types of Monte Carlo moves were used 
for our simulations. The first is a conventional Metropolis 
move, where each particle is displaced in a random direc- 
tion by a random amount. A second Monte Carlo move 
based on the cluster algorithm proposed by Dress and 
Krauth [53] and Liu and Luiijten [2a] was also used. At 
the densities we worked with, neither algorithm is par- 
ticularly efficient in causing very large rearrangements 
in the system configuration. But by mixing two differ- 
ent algorithms that have vastly different properties, we 
hope to minimize equilibration problems characteristic 
of any single algorithm. One Monte Carlo step (MCS) 
in our simulation is defined as having moved each parti- 
cle on the average once using the Metropolis algorithm, 
plus having made one global cluster update. The simula- 
tions reported here were carried out with no fewer than 
5 million MCS for each density. Depending on the equi- 
libration rate, results from the last 1 to 3 million MCS 
were used to collect statistics. Two to four independent 
simulations were carried out for each density for simula- 
tions with a square box, and five to six for those with a 
rectangular box. 

The pressure P was calculated using the virial for- 
mula PAo/NkT = [1 + Trpg(l+)/2]VSp/2, where g(l+) 
is the contact value of the pair correlation function and 
Aq = y/3N/2 is the closed-packed area of the system. 



The calculated pressure P is shown in Fig. 1 as a func- 
tion of density p for N — 512 2 (open circles) and for 
N = 1024 2 (open triangles). Comparing the N = 512 2 
and 1024 2 data to those from Jaster's simulation with 
N = 256 2 , the two sets of data are almost identical for 
p < 0.890, but inside the range p = 0.895 to 0.910, the 
larger size simulations produced a smaller pressure for p 
= 0.895 but larger pressures for p = 0.900 to 0.910. It is 
therefore clear that the apparent van der Waals loop in 
the pressure is a result of finite-size effects, and using even 
larger size simulations, this slight nonmonotonic decrease 
in the pressure should eventually vanish altogether. For 
the N = 1024 2 simulations, the van der Waals loops has 
completely disappeared between p = 0.895 and 0.905, 
with a slight dip in P still visible for p = 0.910. As ex- 
pected, finite-size effects are indeed very pronounced in 
the transition region even for simulations of this magni- 
tude. 

Even though the evidence in Fig. 1 is compelling that 
the van der Waals loop is an artifact of finite-size effects, 
these data alone cannot definitively rule out a first-order 
melting transition. For this, wee need to carefully ana- 
lyze the finite size effects. To disentangle the finite-size 
effects, a detailed scaling analysis must be performed on 
the simulation data. We have found a subblock scaling 
analysis 0, |2(| to be useful for this purpose. With this 
method, a single large size simulation provides informa- 
tion on multiple length scales simultaneously. The sub- 
block scaling analysis was applied to the bond orientation 
order as well as the translational order. 

The bond orientation order is given by i\>\ = 
|(6-/V) _1 ^ z J^j ex P(6*#;j)| 2 , where the sum goes over 
each particles I and its nearest neighbors j and 9ij is the 
angle between the line from I to j with some fixed refer- 
ence axis. According to HNY theory, tpe should have only 
short-range order in the isotropic phase and quasi-long- 
range order in the hexatic phase with exponent rj e < 1/4. 
We calculated tpQ for subblock sizes of Lb = L/64, L/32 
. . . L, where L is the full length of the box and plot the 
results in Fig. 2. For p < 0.895, ip§ clearly scales to zero, 
but for p > 0.900, ip6 appears to scale to a finite value. 

To establish the precise scaling behavior, we plot \nipl 
vs. the natural log of the length of the subblock Lb 
in Fig. 3 for the N = 512 2 simulations. According to 
HNY theory, this plot should show a slope — ry 6 in the 
hexatic phase and —2 in the isotropic phase where there 
is only short-range order. Figure 3 shows that for both 
p = 0.880 (solid triangles) and 0.890 (open diamonds), 
the bond orientation order has no long-ranged correla- 
tions in the long length scale limit, and the size of the 
simulations was large enough to reach the scaling limit. 
We can safely conclude that densities p < 0.890 are in 
the isotropic phase. On the other hand, for the high- 
est densities p — 0.905 (solid diamonds) and 0.910 (open 
squares), the bond orientation order shows an algebraic 
decay with an exponent t/q much smaller than 1/4. This 
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FIG. 2: The bond orientation order parameter ip§ derived 
from the subblock analysis as a function of density for the 
N = 512 2 simulations. For p < 0.895, ipa scales to with 
larger system sizes. For p > 0.905, tpe appears to scale to a 
nonzero value. 
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FIG. 3: Subblock scaling analysis for the bond orientation 
order parameter for the N = 512 2 simulations. The dotted 
line corresponds to a slope of —2 and the dashed line a slope 
of —1/4. The inset shows an expanded view for p — 0.895, 
0.900 and 0.905 in the large length scale region. 

is consistent with the interpretation that these densities 
are either inside the hexatic or the solid phase. 

For the two densities p = 0.895 and 0.900, the interpre- 
tation of the subblock scaling plots is more involved. The 
inset in Fig. 3 shows an expanded view of their behav- 
iors in the large length scale limit. For p = 0.900 (open 
circles) , the bond orientation order shows a slope that is 
very close to —1/4 at large length scales. In the HNY 
scenario, this is consistent with a density inside the hex- 
atic phase, very close to the hexatic-isotropic boundary 
Pi. These evidence suggest that pi<0.900. 

For p = 0.895 (closed squares), the subblock scaling 



plot changes slope twice, first at L/2 and then more 
gradually between L/4 and L/8. The first abrupt slope 
change at L /2 is an artifact of the subblock scaling analy- 
sis which has been discussed by Weber, Marx and Binder 
fll| . The reason for this sudden slope change is that the 
subblocks and the full box actually belongs to two differ- 
ent ensembles - the canonical for the full box and some- 
thing resembling the grand canonical for the subblocks. 
It is therefore possible for the full box to exhibit a differ- 
ent scaling behavior compared to the subblocks when the 
correlation length exceeds the size of the simulation box, 
in which case the full box data point must be excluded 
from the scaling analysis. When this is done, the scaling 
behavior suggests that the orientation order decays alge- 
braically with an exponent larger than 1/4. But clearly 
the scaling limit has not been reached, so it is possible 
that this exponent will continue to increase with lengths 
beyond the size of the present simulation. These evidence 
suggest that p = 0.895 must still be inside the isotropic 
phase but is very close to the isotropic-hexatic boundary. 

Taken together, the pressure data and the subblock 
scaling analysis of the bond orientation order reveal a 
consistent picture. For densities p < 0.895, the system 
is in the isotropic phase. The van der Waals loop in the 
pressure between p = 0.895 and 0.910 observed in previ- 
ous simulations is most certainly due to finite-size effects. 
The bond orientation correlation length increases when 
the isotropic-hexatic boundary pi is approached from be- 
low and it changes from short-range correlation to an al- 
gebraic decay with t]q close to 1/4 at /9,;<0.900, which is 
consistent with previous estimates |l7J. Above pi, the 
exponent t? 6 decreases quickly from 1/4 to zero when 
the hexatic-solid boundary p m is approached from below. 
These findings are consistent with the HNY scenario. 

The fact that 7/ 6 — * for p — > 0.910 has been used 
previously to suggest that the hexatic-solid boundary is 
at p m w 0.910 jllL lri|. The recent data of Jaster, how- 
ever, have placed p m at a much higher value w 0.933 [24| . 
To more accurately locate the hexatic-solid boundary p m , 
we turn to a subblock scaling analysis of the translational 
order i/jf = |iV _1 exp(ik ■ fi)\ 2 , where the wavevector 
k has magnitude 27r/(v / 3/2 ( o) 1 / 2 . In the solid phase, ip* 
is expected to decay algebraically with exponent r\ t = 
1/3. The results for three densities, p = 0.900, 0.910 
and 0.920, are shown in Fig. 4 for the N = 512 2 sim- 
ulations in both a square and a rectangular box with a 
V3 : 2 aspect ratio. For p = 0.900 (triangles) and 0.910 
(squares), the results are consistent with no long-range 
translational order in the large length scale limit for both 
box geometries. This indicates that both of these densi- 
ties are inside the hexatic phase. On the other hand, for 
p = 0.920, the translational order shows apparently dif- 
ferent scaling behaviors for the two box geometries - no 
long-range order in the square box but quasi-long-range 
translational order with an apparent exponent r\ t > 1/3 
for the rectangular box. In fact, comparing the two dif- 
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FIG. 4: Subblock scaling analysis for the translational order 
parameter for the N = 512 2 simulations in a square box and 
a rectangular box. The dotted line corresponds to a slope of 
—2 and the dashed line a slope of —1/3. 



ferent box geometries, we found that the rectangular box 
simulations at this density were much slower to equili- 
brate, leading to the larger error bars on the right panel 
of Fig. 4 for p = 0.920. Since the translational correlation 
length is expected to diverge accordingto HNY theory as 
p approaches the melting density p m |3| , the rectangular 
box used for the simulations at p — 0.920 was probably 
too small to reach the scaling limit. Therefore, we be- 
lieve that p = 0.920 is most likely still inside the hexatic 
phase and has not yet reached the hexatic-solid bound- 
ary. This establishes a lower bound for p m , one that is 
much higher than the value previously suggested [llllr^|. 
But this new lower bound is consistent with the recent 
estimate provided by Jaster based on simulations with 
N up to 1024 2 24]. Since the pressure at p = 0.920 (see 
Fig. 1) is much higher than the pressure inside the appar- 
ent van der Waals loop, this new lower bound for p m also 
provides evidence corroborating the conclusion we have 
drawn from the size-dependence of the pressure-density 
data in Fig. 1 that the apparent van der Waals loop in 
the pressure is not related to a first-order transition. 

In conclusion, we have shown using large-scale com- 
puter simulations with more than a million particles 
that the apparent van der Waals loop observed previ- 
ously in smaller simulations is an artifact of finite-size 
effects. In conjunction with a detailed scaling analysis, 
the data provide compelling evidence for a continuous 
isotropic-hexatic transition as predicted by HNY theory 
at pi<0.900. Scaling analysis of the translational or- 
der also yields a lower bound for the melting density, 
p m > 0.920, one that is much higher than previously 
thought, providing additional evidence that the apparent 



van der Waals loop is not due to a first-order transition. 
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